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ABSTRACT 

Expansions in Chebyshev polynomials are used to study the linear 
stability of one-dimensional magnetohydrodynamic (MHD) quasi-equilibria in the 
presence of finite resistivity and viscosity. The method is modeled on the 
one used by Orszag in accurate computation of solutions of the Orr-Sommerfeld 
equation. Two Reynolds-like numbers involving Alfv&n speeds, length scales, 
kinematic viscosity, and magnetic diffusivity govern the stability boundaries, 
which are determined by the geometric mean of the two Reynolds-like numbers. 
Marginal stability curves, growth rates versus Reynolds-like numbers, and 
growth rates versus parallel wave numbers are exhibited. A numerical result 
which appears general in that instability has been found to be associated with 
inflection points in the current profile, though no general analytical proof 
has emerged. It is possible that nonlinear subcritical three-dimensional 
instabilities may exist, similar to those in Poiseuille and Couette flow. 
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The linear stability of plane shear flows has been one of the most 

intensively studied hydrodynamic problems from the time of Rayleigh, since it 

was thought to hold clues to the nature of turbulence (see, e.g., Lin"'' or 
2 

Maslowe ). Although the linear theory alone appears to be inadequate to predict 

the onset of shear flow instabilities, it remains an important first step in 

any discussion of the problem. We report here on an analogous problem in 

incompressible magnetohydrodynamics (MHD). We report numerical solutions of 

the quiescent-MHD analogue of the Orr-Sommerfeld equation, using spectral 

3 

methods developed by Orszag . 

We begin with the incompressible MHD equations in a familiar dimen- 
sionless form: 

|| = Vx(v x B) + i V 2 B , (1) 

= -y*Vv + B * VJ3 - Vp + V^v , (2) 

supplemented by the conditions that V-v = 0 and V*B = 0. B is the magnetic 
field measured in units of a mean magnetic field magnitude B, say. The velocity 
field is measured in units of the mean Alfven speed = B(l4irp) , where p 
is the mass density, assumed uniform. The dimensionless pressure is p, and it 
is determined by solving the Poisson equation which results from taking the 
divergence of Eq. (2) and using V*8v/3t = 0. The dimensionless numbers S and 
M have the structure of Reynolds numbers. S = C^L/n is the Lundquist number, 
where n is the magnetic diffusivity and L is a macroscopic length scale; 

M = C a L/v is a viscous analogue, where v is a kinematic viscosity. Both ri 
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and v are assumed to be scalars. The regime of most interest is that in 
which S and M are both substantially greater than unity. 

The boundary conditions are taken to be those appropriate to a 
perfectly-conducting, mechanically-impenetrable wall bounding a viscous, 
resistive magnetofluid: v = 0, fi • B = 0, and fi x (V x B) =0, where fi is 

the unit normal at the wall. 

We study the linear stability of the quasi-equilibrium = (B o (y),0,0) 

and = (0,0,0) between parallel, plane infinite boundaries at y = 1 and 

y = -1. The current density is in the z-direction only: j = -DB o , where 

D = d/dy. The configuration described is not a true equilibrium, and the mag- 

netic field will resistively decay according to B Q (y,t) = exp(S tV )B Q (y,0). 

The temporal variation will be assumed to be slow enough to be negligible: 

B Q (y,t) = B Q (y,0) = B Q (y). This implies that our stability boundaries will 

not be accurate in regions of small S; there is in this feature a conceptual 
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difference from the already much studied * * problem of a mean flow parallel 
to a uniform magnetic field with no current, which is a true equilibrium, and 
from Hartmann flow^. 

A linear expansion B = B^°^ + B^ , v = v^ , is assumed, with pro- 
ducts of v^ and B^^ being discarded everywhere in the equations of motion. 
Manipulating the components of the resulting linear equations , we may prove 
a Squire's theorem' 1 ', which implies that for the location of the most unstable 
modes it suffices to consider the two-dimensional case: 3/3z may be set equal 

to zero throughout. All variations with the parallel coordinate x and the 
time t are assumed to be contained in a factor exp( iax-icot) , with a an arbitrary, 
real, parallel wave number and &) = + iw^ a complex eigenvalue. Dahlburg 

and Montgomery have given the eigenvalue equations in the form used here: 
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(D 2 -a 2 ) 2 v = -i(i)M(D 2 -ct 2 )v-iaMB (D 2 -ct 2 )b + iaM(D 2 B )b 

o o 


and 


(D 2 -a 2 + iaiS)b = -ictSB v. 


(M 


Here b and v are the y components of B^ and v^^, and depend only upon y. 

The boundary conditions become V = 0, Dv = 0, and b = 0 at y = 1 and y = -1. 

Equations (3) and (U) are the magnetostatic analogue of the Orr- 

1 2 3 2 2 2 

Sommerfeld Equation, which in the same notation * ’ is (D -a ) v = 

2 2 2 

iaR[(U o -w/a) (D -a )v-(D U q )v], where U Q (y) is a shear flow velocity profile in 
the x direction, R is the Reynolds number, and the boundary conditions are 
that v = 0, Dv = 0 at y = ±1. 

Equations (3) and (k) are quite similar to eigenvalue problems 
arising in connection with confinement of thermonuclear plasmas. The litera- 
ture of "tearing modes" is extensive, and we may cite the central papers of 

Furth, Killeen, and Rosenbluth®, of Wesson^, of Coppi, Greene, and Johnson^, 
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of Furth, Rutherford, and Selberg , and of Dibiase and Killeen . Concern 

has frequently been with the non-viscous (M = °°) case, which lowers the order 

of the differential equations. Viscous results from a linear initial-value 

12 

computation have been reported by Dibiase and Killeen for the compressible 
case, and to the extent that the results can be compared, ours do not appear 
to disagree with theirs. Because for plasmas of interest to date, the calcu- 
lated viscosity coefficients give estimates of V at least as great as those 

13 

for r) (see, e.g. , Braginskii ), it seems desirable to retain viscous effects 

even at the price of reusing the order of the eigenvalue problem to the point 

where results can be extracted only numerically. 

Analytical information is difficult to extract even for the simpler 

2 1 ^ 

Orr-Sommerfeld equation (e.g., Maslowe or Reid ) and numerical solution is 
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indicated here. We use a spectral technique closely patterned on the method 
3 

used by Orszag to calculate critical Reynolds numbers for Poiseuille flow to 
six-figure accuracy. It is to be expected that spectral methods will find 
further applications in plasma physics beyond the immediate ones. 

The mean magnetic field B Q (y) and the perturbation quantities v and 
b are expanded in truncated Chebyshev series 

N 

B o<*> - E ®„ T „(y) 

n=0 

N 

v(y) = Z v T n (y) 
n=0 

N 

B(y) = s b T (y), (5) 

n=0 n n 

where T Q (y) is the n Chebyshev polynomial of the first kind, and B r , v^, 
b are the respective expansion coefficients. 

The equations satisfied by the (unknown) expansion coefficients are 
obtained by substituting the N -*■ 00 expansions of (5) into Eqs. (3) and (4), 
each of which produces a countably infinite number of equations in the expan- 
sion coefficients for n = 0, 1, 2, ...» when the orthogonality and recursion 

3 

relations are used. We then set all coefficients beyond n = N to zero, lose 

the n = 0 to n = N-4 equations from (3), the n = 0 to N-2 equations from (4), 

N N N _ 

and the boundary conditions E v = 0, E (-1) v =0, E n v =0, 

N N n=0 n u n=0 n n=l n 

E (-l) n n^v =0, E b =0, and E (-l) n b =0. This method of truncation 
n=l n=0 n=0 

is called the "tau approximation" after Lanczos ^ ; Gottlieb and Orszag give 
a general discussion of the use of the tau method in the Chebyshev case. 

This spectral discretization process yields a generalized eigenvalue 
problem which can be written els Ax = coBx, where the vector x = (v q , v^, ... v^, 
b Q , b^, ... b^), and A and B axe non-symmetric (2N+2) by (2N+2) square matrices. 
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As is customary for this type of stability problem, either global 
or local, methods are used to determine the eigenvalues. The global method 
is based on the QR algorithm (Wilkinson - , Gary and Helgason xo ) and produces 
a full spectrum of eigenvalues. It is employed when no good guess for the 
least stable (or most unstable) eigenvalue is available. The local method 
employs inverse Rayleigh iteration ' and converges to the eigenvalue (and its 
associated eigenfunction) closest to the initial guess for the eigenvalue. 

The global method may be used to identify the eigenvalue with the largest 
imaginary part, and the local method is useful when making a series of computa- 
tions in which either the wavenumber or the Reynolds numbers are slowly varied. 

For functions B Q (y) which are antisymmetric about y = 0, it is readily 
inferred from Eqs. (3) and (1») that v and b are of opposite parity when re- 
flected about y = 0. We have confined attention to the case B Q (y) = -B^t-y) 
with an associated current distribution J q = -DB q which has even parity about 
y = 0: the classic "sheet pinch" configuration. This configuration (indeed, 

any B (y) profile) can be rigorously proved to be stable in the ideal limit 
(M = °°, S = °°); any instabilities must result from finite values of S or M 
or (in our case) both. 

We have solved for the several eigenfunctions corresponding to the 
largest values of for four different antisymmetric profiles B Q (y): 

B o I (y) = y - y 3 /3 
B^fy) = tan -1 yy - yy(Y 2 +l) -1 
B q (y) = y - y /2i 
B Q IV (y) = sinh -1 yy - yy( Y^l) -1 ^ 2 . 


( 6 ) 
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Two numerical results have characterized all runs performed, and 
though ve have been unable to prove either one analytically, we suspect they 
are generally true: (1) as S or M is raised a first -unstable eigenvalue 
always appears (ok > 0) at finite values of S and M with = 0; and (2) a 
necessary condition that instabilities appear is that the current profile J Q 
shall have an inflection point between y = -1 and y = +1. Steep current gradi- 
ents alone seem insufficient to produce instability. For example, B^ (y), 
which has a large maximum current gradient of 20 near the walls, was found to 
be stable up to M = 10 \ S = 10^ (for a = l) , while profiles such as which 
did contain inflection points would characteristically be unstable for S and M 
no greater than a few tens, with considerably smaller current gradients involved. 

Particularly extensive investigations were carried out for the profile 
B^(y) for various values of a, S, M, and the "stretching parameter" y. Figure 
1 is a plot of B^(y) and its associated current profile j^(y) = -DB^ as a 
function of y for y = 10. This case will be used to illustrate the results in 
Figs. 2 through 9- 

Figure 2 shows typical eigenfunctions , stable and unstable, computed 
from the B* 1 of Fig. 1. Figure 2a applies to S = M = 10, for the least-damped 
eigenvalue on = -0.1695- For this case, u=0, o=1.0, b= ib^ is purely 
imaginary, and v = v^ 'is purely real. This last property always applies to the 
eigenfunction with the greatest OK. Figure 2b shows b^ , v^ for an eigenfunc- 
tion immediately above the instability threshold, with S = M = 20, 
ok = 0. 069 1+0, a = 1.0, = 0. Figure 2c shows K , v r well above the thres- 

hold, with S = M = 1000, a = 1.0, ok = 0.19687, and o> r - 0. Figure 2d shows 
the case S = 10, M = 10^, with a = 1.0, ok = 0.1+397 and. a> r = 0. Figure 2e 
shows the case S = 10^, M = 10, GK = 0. 002537 > ot = 1.0, = 0, and 
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illustrates the (perhaps unsurprising) result that viscosity is better at 
suppressing unstable growth than resistivity. 

The damped modes have, in general, both 63^ and 03^ finite. Typical 
eigenfunctions for a highly damped mode for are shown in Fig. 3 
(o) r = -0.U0U92, ok = -6.01303, S = M = 20, a = 1.0). 

- 1/2 

At the stability threshold 63 = 0, the scaling v' = vM , 
b' = bS~ ' reduces Eqs. (3) and (1+) to a pair of equations for v' and b' 
which depend upon S and M only in the combination v^SM (and, of course, upon 
a). The neutral stability curve in the SM plane, across which 63^ becomes 
positive for some a, is therefore a hyperbola, approximately SM = 231.9* 

The computed location of this hyperbola (for with y = 10) is shown in 
Fig. 4. The relatively low values of the critical Reynolds numbers (two orders 
of magnitude or more below the corresponding hydrodynamic ones for shear flows) 
are perhaps the most significant feature of this graph. It is also interest- 
ing that for low enough values of either Reynolds number, stability will always 

result for any fixed value of the other, but since SM increases with tempera- 

13 

ture, according to kinetic theory estimates , at high enough temperatures, 

we may always expect to be on the unstable side of the boundary. 

To make a comparison with traditional - * - hydrodynamic plots. Fig. 5 

exhibits a set of marginal stability curves: o> = 0 in the a, S plane for 

fixed M for B^ 1 and y = 10. The first unstable a (which, for the reasons 

previously noted, must be the same for all such curves) is a = = l.l8H±0.005. 

We have generated the same curves (not shown) in the a, M plane for fixed S. 

Figure 6a shows a growth rate (o)^ vs. S) curve for M = 10. Figure 

-3/5 

6b shows a logarithmic plot for large values of S, illustrating the S 

g 

regime identified analytically by Furth et al . Figure 7 shows the somewhat 
different behavior of 63^ vs. M for fixed S. Figure 8a shows a contour 
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plot of two periods of the magnetic field lines for the B** equi- 

librium plus a 20 percent admixture of the eigenfuction shown in Fig. 2b. Figure 
8b shows the streamlines of the velocity field for the same eigenfuction. 

Figure 9 shows the growth rate ok vs', a for several values of S and M. 

The results presented have all been computational and are not the 
result of an asymptotic "tearing layer" analysis. The marginal stability 
curves such as Fig. 5 or the stability hyperbola of Fig. U could only have 
been obtained numerically. Despite these new results, we are well aware of 
other potentially important physical effects which have been left out of the 
analysis, such as compressibility, or the effects of spatially varying vis- 
cosity and resistivity tensors. Nevertheless, we believe there to be some 
value in bringing a more classical hydrodynamic perspective to bear on this 
highly idealized problem. It is not impossible that subcritical instability 
thresholds from nonlinear three-dimensional computations may be identified 
for this problem an they have been for shear flows (see, e.g. , Orszag and 
Kells'^ or Orszag and Patera^). Further computations of a more elaborate 

kind will be required to test this possibility; for further discussion of the 

21 

hydrodynamic antecedents, see the review by Herbert 
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FIGURE CAPTIONS 

Fig. 1 B q = B^(y) and its associated current profile = -DB q , for y = 10. 

Fig. 2a Eigenfunctions b = iK , v = v r for the least-damped eigenvalues for 

M = 10, S = 10. 


Fig. 2b 

Fig. 2c 
Fig. 2d 
Fig. 2e 
Fig. 3 
Fig, b 

Fig. 1 
Fig. 6 a 


Eigenfunctions b = ib^ , v = v r slightly above the instability thres- 
hold: M = 20, S = 20. 

Unstable eigenfunctions b = ib^ , v = v^ for S = M = 1000. 

Unstable eigenfunctions b = iIk , v = v r for S = 10, M = 10^. 

Unstable eigenfunctions b = itK , v = v^ for S = 10^, M = 10. 

Highly damped complex eigenfunctions for M = S = 20. 

Locus of critical Reynolds-like numbers S = S c , M = in the M3 
plane, as determined by computation. 

The neutral stability curve u = 0 in the a, S plane for M = 1, 10, 
and 1000 for B^ 1 and y = 10. 

Growth rates (gk vs. S) for a = 1.0, (on B* 1 with y = 10), 

M = 10, 100, and 1000. 


Fig. 6b Growth rate ok vs. S for fixed M, a = 1.0, B q with y = 10. 

Fig. 7 Growth rate on vs. M for fixed S = 10, 100, 1000; a = 1.0, y = 10 
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Fig. 8a Contour plot of magnetic field lines for a typical unstable eigen- 
function near the threshold (S = M = 20). Field lines are equilibrium 
plus 20 percent admixture of eigenfunction. 

Fig. 8b Velocity streamlines corresponding to eigenfunction represented in 
Fig. 8a. 

Fig. 9 Growth rate vs. parallel wave number a for , M = S = 100; 

M = S = 500; M = 1000, S = 100; M = 100, S = 1000. 
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Fig. 1 B = B (y) and its associated current profile J = -DB , for y = 10. 
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Fig. 3 Highly damped complex eigenfunctions for 
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Fl,^. 8a Contour plot of magnetic field lines for a typical unstable eigen- 

function near the threshold (S = M = 20 ). Field lines are. equilibrium 
plus -20 percent admixture of eigenfunction. 
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FiR. 8b Velocity streamlines corresponding to eigenfunction represented in 


Fig. 8a. 
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